library(survey)
library(foreign)
library(ggplot2)
library(cowplot)
library(dplyr)
library(egg)
library(ggpubr)
library(gridGraphics)

## see wave1_CLEANING.R for data cleaning
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave1_data/clean_dataw1.RData")

## Create weighted survey design object
w1_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

addline_format <- function(x,...){
  gsub('\\s','\n',x)
}

#### TREATMENT EFFECTS FOR CONFIDENCE IN MILITARY ####

treax_df <- data.frame(otva = c(svyttest(Q8_d ~ ASSIGNMENT_A_2, w1_design)$estimate,
                                svyttest(Q8_d ~ ASSIGNMENT_A_3, w1_design)$estimate,
                                svyttest(Q8_d ~ ASSIGNMENT_A_4, w1_design)$estimate),
                       ci_lo = c(svyttest(Q8_d ~ ASSIGNMENT_A_2, w1_design)$conf.int[1],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_3, w1_design)$conf.int[1],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_4, w1_design)$conf.int[1]),
                       ci_hi = c(svyttest(Q8_d ~ ASSIGNMENT_A_2, w1_design)$conf.int[2],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_3, w1_design)$conf.int[2],
                                 svyttest(Q8_d ~ ASSIGNMENT_A_4, w1_design)$conf.int[2]),
                       cond = c(1:3))

fig71_a <- ggplot(treax_df) + 
  geom_pointrange(aes(x = cond, y = otva, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = otva, label = sprintf("%0.2f", round(otva, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1:3),
                     labels = addline_format(c("Republican Military", "Democrat Military",
                                               "Polarized Military"))) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  geom_hline(yintercept = 0, colour = "red", lty = 1) + theme_bw()
fig71_a
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_diff_treateffects.jpeg", 
       plot = fig71_a, dpi = 1000)

#### TREATMENT EFFECTS BY CONFIDENCE BY PARTY ####

treax2_df <- data.frame(rep = c(svyttest(Q8_d ~ A2_rep, w1_design)$estimate,
                                svyttest(Q8_d ~ A3_rep, w1_design)$estimate,
                                svyttest(Q8_d ~ A4_rep, w1_design)$estimate),
                        dem = c(svyttest(Q8_d ~ A2_dem, w1_design)$estimate,
                                svyttest(Q8_d ~ A3_dem, w1_design)$estimate,
                                svyttest(Q8_d ~ A4_dem, w1_design)$estimate),
                        ind = c(svyttest(Q8_d ~ A2_ind, w1_design)$estimate,
                                svyttest(Q8_d ~ A3_ind, w1_design)$estimate,
                                svyttest(Q8_d ~ A4_ind, w1_design)$estimate),
                        rep_lo = c(svyttest(Q8_d ~ A2_rep, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A3_rep, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A4_rep, w1_design)$conf.int[1]),
                        rep_hi = c(svyttest(Q8_d ~ A2_rep, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A3_rep, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A4_rep, w1_design)$conf.int[2]),
                        dem_lo = c(svyttest(Q8_d ~ A2_dem, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A3_dem, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A4_dem, w1_design)$conf.int[1]),
                        dem_hi = c(svyttest(Q8_d ~ A2_dem, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A3_dem, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A4_dem, w1_design)$conf.int[2]),
                        ind_lo = c(svyttest(Q8_d ~ A2_ind, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A3_ind, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A4_ind, w1_design)$conf.int[1]),
                        ind_hi = c(svyttest(Q8_d ~ A2_ind, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A3_ind, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A4_ind, w1_design)$conf.int[2]),
                        dcond = c(1,5,9),
                        icond = c(2,6,10),
                        rcond = c(3,7,11))

fig71_b <- ggplot(treax2_df) + 
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10),
                     labels = addline_format(c("Republican Military", "Democrat Military",
                                               "Polarized Military"))) +
  ylab("Difference in Means from Control Condition") +
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  theme_bw()
fig71_b
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_diff_treatparty.jpeg", 
       plot = fig71_b, dpi = 1000)

#### TREATMENT EFFECTS - NEGATIVE PARTISANSHIP ####

copart1 <- svyby(~Q8_d, ~A_copartisan, w1_design, svymean, na.rm = TRUE)

np_df1 <- data.frame(outs = c(NA, copart1[,2], NA),
                     ci_lo = c(NA, copart1[,2] - 1.96*copart1[,3], NA),
                     ci_hi = c(NA, copart1[,2] + 1.96*copart1[,3], NA),
                     cond = c(1.5,2,3,4,4.5))

fig7_np1 <- ggplot(np_df1) + 
  geom_pointrange(aes(x = cond, y = outs, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = outs, label = sprintf("%0.2f", round(outs, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("") + 
  scale_x_continuous(breaks = c(2:4),
                     labels = addline_format(c("Control", "Co-Partisan", "Cross-Partisan"))) +
  ylab("Respondent Confidence in Military") + ylim(0.4,0.90) +
  ggtitle("Wave 1") +
  theme_bw()
fig7_np1
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_np.jpeg", 
       plot = fig7_np1, dpi = 1000)

#### TREATMENT EFFECTS - NEGATIVE PARTISANSHIP, DIFFERENCE-IN-MEANS ####

copart1_df <- data.frame(otva = c(NA, svyttest(Q8_d ~ A_copartisanc, w1_design)$estimate,
                                  svyttest(Q8_d ~ A_copartisanx, w1_design)$estimate, NA),
                         ci_lo = c(NA, svyttest(Q8_d ~ A_copartisanc, w1_design)$conf.int[1],
                                   svyttest(Q8_d ~ A_copartisanx, w1_design)$conf.int[1], NA),
                         ci_hi = c(NA, svyttest(Q8_d ~ A_copartisanc, w1_design)$conf.int[2],
                                   svyttest(Q8_d ~ A_copartisanx, w1_design)$conf.int[2], NA),
                         cond = c(0.5,1,2,2.5))

fignp_d1 <- ggplot(copart1_df) + 
  geom_hline(yintercept = 0, colour = "red", lty = 1) + 
  geom_pointrange(aes(x = cond, y = otva, ymin = ci_lo, ymax = ci_hi), fatten = 2, size = 0.5) +
  geom_text(aes(x = cond, y = otva, label = sprintf("%0.2f", round(otva, digits = 2))), 
            size = 3,
            vjust = 2.5) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1,2),
                     labels = c("Co-Partisan", "Cross-Partisan")) +
  ylab("Difference in Means from Control Condition") + ylim(-0.2, 0.2) +
  ggtitle("Wave 1") + theme_bw()
fignp_d1
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_np_diff.jpeg", 
       plot = fignp_d1, dpi = 700)


#### TREATMENT EFFECTS FOR BORDER DEPLOYMENT BY PARTY ####

treat3_df <- data.frame(rep = c(svyttest(Q42A_d ~ C3_rep, w1_design)$estimate,
                                svyttest(Q42A_d ~ C4_rep, w1_design)$estimate,
                                svyttest(Q42A_d ~ C5_rep, w1_design)$estimate,
                                svyttest(Q42A_d ~ C6_rep, w1_design)$estimate),
                        dem = c(svyttest(Q42A_d ~ C3_dem, w1_design)$estimate,
                                svyttest(Q42A_d ~ C4_dem, w1_design)$estimate,
                                svyttest(Q42A_d ~ C5_dem, w1_design)$estimate,
                                svyttest(Q42A_d ~ C6_dem, w1_design)$estimate),
                        ind = c(svyttest(Q42A_d ~ C3_ind, w1_design)$estimate,
                                svyttest(Q42A_d ~ C4_ind, w1_design)$estimate,
                                svyttest(Q42A_d ~ C5_ind, w1_design)$estimate,
                                svyttest(Q42A_d ~ C6_ind, w1_design)$estimate),
                        rep_lo = c(svyttest(Q42A_d ~ C3_rep, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C4_rep, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C5_rep, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C6_rep, w1_design)$conf.int[1]),
                        rep_hi = c(svyttest(Q42A_d ~ C3_rep, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C4_rep, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C5_rep, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C6_rep, w1_design)$conf.int[2]),
                        dem_lo = c(svyttest(Q42A_d ~ C3_dem, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C4_dem, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C5_dem, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C6_dem, w1_design)$conf.int[1]),
                        dem_hi = c(svyttest(Q42A_d ~ C3_dem, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C4_dem, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C5_dem, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C6_dem, w1_design)$conf.int[2]),
                        ind_lo = c(svyttest(Q42A_d ~ C3_ind, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C4_ind, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C5_ind, w1_design)$conf.int[1],
                                   svyttest(Q42A_d ~ C6_ind, w1_design)$conf.int[1]),
                        ind_hi = c(svyttest(Q42A_d ~ C3_ind, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C4_ind, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C5_ind, w1_design)$conf.int[2],
                                   svyttest(Q42A_d ~ C6_ind, w1_design)$conf.int[2]),
                        dcond = c(1,5,9,13),
                        icond = c(2,6,10,14),
                        rcond = c(3,7,11,15))

fig71_c <- ggplot(treat3_df) + 
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10, 14),
                     labels = addline_format(c("Trump Politicizes", 
                                               "Democrat Politicizes",
                                               "Military Supports",
                                               "Military Opposes"))) +
  ylab("Difference in Means from Control Condition") +
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  theme_bw()
fig71_c
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_diff_border.jpeg", 
       plot = fig71_c, dpi = 1000)

#### TRUST BY BORDER DEPLOYMENT AND PARTY ####

treat4_df <- data.frame(rep = c(svyttest(Q42F_d ~ C3_rep, w1_design)$estimate,
                                svyttest(Q42F_d ~ C4_rep, w1_design)$estimate,
                                svyttest(Q42F_d ~ C5_rep, w1_design)$estimate,
                                svyttest(Q42F_d ~ C6_rep, w1_design)$estimate),
                        dem = c(svyttest(Q42F_d ~ C3_dem, w1_design)$estimate,
                                svyttest(Q42F_d ~ C4_dem, w1_design)$estimate,
                                svyttest(Q42F_d ~ C5_dem, w1_design)$estimate,
                                svyttest(Q42F_d ~ C6_dem, w1_design)$estimate),
                        ind = c(svyttest(Q42F_d ~ C3_ind, w1_design)$estimate,
                                svyttest(Q42F_d ~ C4_ind, w1_design)$estimate,
                                svyttest(Q42F_d ~ C5_ind, w1_design)$estimate,
                                svyttest(Q42F_d ~ C6_ind, w1_design)$estimate),
                        rep_lo = c(svyttest(Q42F_d ~ C3_rep, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C4_rep, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C5_rep, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C6_rep, w1_design)$conf.int[1]),
                        rep_hi = c(svyttest(Q42F_d ~ C3_rep, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C4_rep, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C5_rep, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C6_rep, w1_design)$conf.int[2]),
                        dem_lo = c(svyttest(Q42F_d ~ C3_dem, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C4_dem, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C5_dem, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C6_dem, w1_design)$conf.int[1]),
                        dem_hi = c(svyttest(Q42F_d ~ C3_dem, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C4_dem, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C5_dem, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C6_dem, w1_design)$conf.int[2]),
                        ind_lo = c(svyttest(Q42F_d ~ C3_ind, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C4_ind, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C5_ind, w1_design)$conf.int[1],
                                   svyttest(Q42F_d ~ C6_ind, w1_design)$conf.int[1]),
                        ind_hi = c(svyttest(Q42F_d ~ C3_ind, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C4_ind, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C5_ind, w1_design)$conf.int[2],
                                   svyttest(Q42F_d ~ C6_ind, w1_design)$conf.int[2]),
                        dcond = c(1,5,9,13),
                        icond = c(2,6,10,14),
                        rcond = c(3,7,11,15))

fig71_d <- ggplot(treat4_df) + 
  geom_pointrange(aes(x = rcond, y = rep, ymin = rep_lo, ymax = rep_hi), color = "red") +
  geom_pointrange(aes(x = icond, y = ind, ymin = ind_lo, ymax = ind_hi), color = "gray60") +
  geom_pointrange(aes(x = dcond, y = dem, ymin = dem_lo, ymax = dem_hi), color = "blue") +
  geom_text(aes(x = rcond, y = rep, label = sprintf("%0.2f", round(rep, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = icond, y = ind, label = sprintf("%0.2f", round(ind, digits = 2))), size = 3,
            vjust = 2) +
  geom_text(aes(x = dcond, y = dem, label = sprintf("%0.2f", round(dem, digits = 2))), size = 3,
            vjust = 2) + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10, 14),
                     labels = addline_format(c("Trump Politicizes", 
                                               "Democrat Politicizes",
                                               "Military Supports",
                                               "Military Opposes"))) +
  ylab("Difference in Means from Control Condition") +
  geom_hline(yintercept = 0, colour = "black", lty = 1) +
  theme_bw()
fig71_d
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 6 Figures/wave1_diff_bordertrust.jpeg", 
       plot = fig71_d, dpi = 1000)
